Using only tracks that passed QC


Figure 1 - Baseline Dominance

Baseline Cumulative DS

  • Cumulative David’s Score based on chases

Correlations with SB behaviors

  • Cumulative David’s Scores from days 1:4
  • SB behaviors summarized as means of days 2:4

Stats

Social box features with significant q-values:

Female Male
MedianContactDuration 0 0
MeanContactDuration 0 0
FractionOfTimeInTheOpenOutside 0 0
MeanSpeedOutside 0 0
FractionOfTimeInSmallNestOutside 0 0
FractionOfTimeInFeederOutside 0 0
FoodOrWaterPerTimeOutside 0 0
MedianSpeedOutside 0 0
FractionOfTimeOnRampOutside 0 0
DistanceFromWallsInOpen 0 0
GridMI6x6 0 0
EntropyOutside 0 0
DistanceFromNest 0 0
HighPlacePerTimeOutside 0 0
ContactRate 0 0
VisitsOutsideRate 0 0
BeingApproachedRate 0 0
ForagingCorrelation 0 0
AngularVelocity 0 0
FractionOfTimeAtHighPlace 0 0
FractionOfTimeInWaterOutside 0 1
ApproachRateOutside 1 0
TangentialVelocity 0 0
DiffBetweenApproachesAndChases 0 0
FractionOfTimeInContactOutside 0 0
ContactRateOutside 0 1
FractionOfTimeNearFoodOrWater 0 1
DistanceOutside 0 1
FractionOfTimeInLabyrinthOutside 1 0
NumberOfApproachs 0 1
ApproachRate 0 1
FractionOfTimeAloneOutside 0 1
ProximateVsDistantFood 0 1
GridEntropy6x6 0 1
Entropy 0 1
NumberOfApproachesPerMiceOut 0 1
FractionOfTimeOutside 0 1
BeingApproachedRateOutside 0 1
FractionOfBeingApproachedPerContact 0 1
FractionOfApproachEscapeBehaviorPerAggression 0 1
AggressiveEscapeRate 0 1
FractionOfApproachesPerContact 0 1
BeingFollowedRate 1 1
NAEscapeRate 1 1
AggressiveChaseRateOutside 1 1
AggressiveEscapeRateOutside 0 1
FractionOfEscapesPerContact 0 1
AggressiveChaseRate 1 1
NAEscapeRateOutside 1 1
FractionOfChasesPerContact 1 1
NAChaseRateOutside 1 1
FollowRateOutside 1 1
BeingFollowedRateOutside 1 1
FollowRate 1 1
NAChaseRate 1 1
FractionOfBeingFollowedPerContact 1 1
FractionOfNAEscapesPerContact 1 1
FractionOfFollowsPerContact 1 1
FractionOfNAChasesPerContact 1 1
ProximateVsDistantWater 0 0

Distance Outside:

Fraction of Time Outside:

Grid Entropy [6x6]

DS time correlations

  • Pearson correlations of DS values between Day 1 and the following days

Sum of daily changes

# Generate a distribution of expected rank changes
set.seed(42)
sampleGroupRanks = function (i) {sample(1:i, size = i, replace = F)}

# Fake data
nMice = length(unique(obj$master$RealMouseNumber))
nDays = 4
nIter = 10

ref = lapply(1:nIter, function(i) {
  dat = data.frame(
    Day = rep(1:nDays, nMice),
    MouseNumber = rep(1:nMice, each = nDays), 
    GroupNumber = rep(1:(nMice/4), each = 4*nDays)
  )
  dat = dat %>% group_by(GroupNumber, Day) %>% 
    nest() %>%
    mutate(
      Rank = map(data, function(d) sampleGroupRanks(nrow(d)))
    ) %>% unnest()
  
  d = dat %>% group_by(GroupNumber, MouseNumber) %>% 
    arrange(GroupNumber, MouseNumber, Day) %>% 
    mutate(
      lagRank = lag(Rank, 1),
      RankChange = abs(Rank - lagRank)
    ) %>% filter(!is.na(lagRank)) 
  ref = d %>% group_by(GroupNumber, MouseNumber, RankChange) %>%
    tally() %>%
    group_by(GroupNumber, MouseNumber) %>% 
    dplyr::mutate(
      s = sum(n),
      SumRankChangesFraction = n / sum(n) 
    ) %>% 
    group_by(RankChange) %>% 
    tally() %>% 
    dplyr::mutate(
      SumRankChangesFraction = n / sum(n)
    )
  ref$iter = i
  return(ref)
})

ref = do.call('rbind', ref)
ref = ref %>% group_by(RankChange) %>% 
  dplyr::summarise(
    M = mean(SumRankChangesFraction, na.rm = T), 
    Qt = qt(0.975,n()-1) * sd(SumRankChangesFraction)/sqrt(n()),
    CI_95 = M + Qt,
    CI_05 = M - Qt
  )

# Real final ranks
r = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(MouseID)) %>% rename(DS_cum = DS)
r = left_join(r, obj$mice[,c('MouseID', 'GroupNumber', 'Sex', 'RealMouseNumber')])

r = r %>% group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-DS_cum), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )

# Real ranks throughout
p = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = T) %>% 
  dplyr::select(GroupNumber, Sex, Condition, Day, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(MouseID)) %>% 
  group_by(GroupNumber, Sex, Condition, Day) %>% 
  mutate(DS_rank = rank(-DS))
p = left_join(p, obj$mice[,c('RealMouseNumber', 'GroupNumber', 'Sex', 'MouseID')])
p$Stage = ifelse(p$Day %in% 1:4, 'Baseline', 'Stress')
p = left_join(p, r[,c('RealMouseNumber', 'Rank')])
p = p %>% group_by(GroupNumber, Sex, Condition, MouseID, Rank) %>% 
  arrange(GroupNumber, Sex, MouseID, Day) %>% 
  mutate(
    lagRank = lag(DS_rank),
    RankChange = abs(DS_rank - lagRank)
  ) %>% filter(!is.na(RankChange))

pp = p %>% 
  # group_by(Sex, Rank, RankChange, RealMouseNumber) %>% 
  group_by(Sex, Rank, RankChange) %>%
  tally() %>%
  # group_by(Sex, Rank, RealMouseNumber) %>% 
  group_by(Sex, Rank) %>%
  mutate(
    s = sum(n),
    SumRankChangesFraction = n / sum(n, na.rm = T) 
  ) 

ggplot(pp, aes(x = RankChange, y = SumRankChangesFraction)) +
  geom_ribbon(data = ref, aes(ymax = CI_95, ymin = CI_05, y = M), 
              fill = 'grey', alpha = 0.5) +
  geom_line(data = ref, color = 'grey', lty = 1, aes(y = M)) +
  geom_line(aes(col = Rank)) +
  geom_point(aes(fill = Rank), pch = 21, color = 'grey20', size = 2.5) +
  facet_grid(.~Sex) + 
  labs(y = 'Probability') +
  scale_color_manual(values = obj$colors_rank) +
  scale_fill_manual(values = obj$colors_rank) + 
  theme(legend.position = 'top')

Change/no change version

# Real final ranks
r = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(MouseID)) %>% rename(DS_cum = DS)
r = left_join(r, obj$mice[,c('MouseID', 'GroupNumber', 'Sex', 'RealMouseNumber')])

r = r %>% group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-DS_cum), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )

# Real ranks throughout
p = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = T) %>% 
  dplyr::select(GroupNumber, Sex, Condition, Day, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(MouseID)) %>% 
  group_by(GroupNumber, Sex, Condition, Day) %>% 
  mutate(DS_rank = rank(-DS))
p = left_join(p, obj$mice[,c('RealMouseNumber', 'GroupNumber', 'Sex', 'MouseID')])
p$Stage = ifelse(p$Day %in% 1:4, 'Baseline', 'Stress')
p = left_join(p, r[,c('RealMouseNumber', 'Rank')])

p = p %>% group_by(GroupNumber, Sex, Condition, MouseID, Rank) %>% 
  arrange(GroupNumber, Sex, MouseID, Day) %>% 
  mutate(
    lagRank = lag(DS_rank),
    RankChange = abs(DS_rank - lagRank)
  ) %>% filter(!is.na(RankChange))

pp = p %>% 
  dplyr::mutate(
    NoChange = RankChange == 0
  )

ppp = pp %>% group_by(Rank, Sex) %>%
  dplyr::summarise(
    N = n(), 
    S = sum(NoChange)
  ) %>%
  dplyr::mutate(
    P = S / N
  )

ppp$sig = NA
ppp$sig[ppp$Rank == 'alpha'] = '***'
ppp$sig[ppp$Rank == 'delta' & ppp$Sex == 'Male'] = '***'


ppp$Rank = plyr::mapvalues(ppp$Rank, c('alpha', 'beta', 'gamma', 'delta'), 
                           c( '\u03B1', '\u03B2', '\u03B3', '\u03B4'))
f1_3 = ggplot(ppp, aes(x = Rank, y = P / .25, fill = Rank)) +
  geom_bar(stat = 'identity', width = 0.7, alpha = 1, aes(color = Rank)) +
  geom_hline(yintercept = 1, lty = 1, color = 'grey30') +
  geom_text(aes(label = Rank, y = 1), nudge_y = 0.08, parse = F, fontface = 'italic', color = 'grey30') +
  geom_text(aes(label = N), color = 'grey20', nudge_y = 0.08, parse = F, size = 3, fontface = 'italic', color = 'grey30') +
  geom_text(aes(label = sig), color = 'grey20', nudge_y = 0.15, parse = F, size = 4, fontface = 'bold', color = 'grey30') +
  facet_grid(.~Sex) + 
  labs(y = 'Rank maintenance odds', x = '') +
  scale_color_manual(values = obj$colors_rank) +
  scale_fill_manual(values = obj$colors_rank) +
  scale_y_continuous(trans = 'log2', breaks = c(0.66, 1, 1.5, 2, 2.5, 3), limits = c(0.66, 3)) +
  theme(legend.position = 'none', 
        axis.text.x =  element_blank())

f1_3

Stats:

DS baseline v stress

Baseline v. Stress DS

Stats

Males:

Females

Chases after stress

Chases from baseline to stress

  • Mean number of chases daily: baseline (days 2:4) vs. stress
pp = h2 %>% group_by(GroupNumber, MouseID, Sex, Condition) %>% 
  summarise_at(.vars = c('wins', 'wins_stress'), mean, na.rm = T) %>% 
  gather(Stress, Chases, -GroupNumber, -Sex, -Condition, -MouseID) %>% 
  mutate(Stress = plyr::mapvalues(Stress, c('wins', 'wins_stress'), c('Baseline', 'Acute Stress')), 
         Stress = factor(Stress, levels = c('Baseline', 'Acute Stress')))

Md = pp %>% group_by(Sex, Stress) %>% 
  dplyr::summarise(
    Md = median(Chases, na.rm = T),
    Mn = mean(Chases, na.rm = T),
    SE = sd(Chases, na.rm = T) / sqrt(n()),
    MAD = stats::mad(Chases, na.rm = T)
  )
# ppp = left_join(pp, Md)
f1_5 = ggplot(pp,
              aes(x = Chases, fill = Sex, color = Sex)) +
  geom_density(alpha = 0.7) +
  geom_rug(aes(color = Sex), size = 1, show.legend = F) +
  geom_errorbarh(data = Md[Md$Sex == 'Male',], inherit.aes = F, 
                 aes(xmin = Md, xmax = Md + MAD, y = -0.1, col = Sex), height = 0.1, show.legend = F) +
  geom_point(data = Md[Md$Sex == 'Male',], 
             pch = 21, size = 2.5, aes(x = Md, y = -0.1), color = 'grey20', show.legend = F) +
  geom_errorbarh(data = Md[Md$Sex == 'Female',], inherit.aes = F, 
                 aes(xmin = Md, xmax = Md + MAD, y = -0.2, col = Sex), height = 0.1, show.legend = F) +
  geom_point(data = Md[Md$Sex == 'Female',], 
             pch = 21, size = 2.5, aes(x = Md, y = -0.2), color = 'grey20', show.legend = F) +
  scale_fill_manual(values = obj$colors_sex) +
  scale_color_manual(values = obj$colors_sex) +
  facet_grid(Stress~.) + 
  scale_x_log10() + labs(y = '') +
  theme(legend.position = c(0.1, 0.9), 
        legend.title = element_blank(),
        panel.grid.major.x = element_line(color = 'grey', size = 0.2, linetype = 2),
        # panel.grid.minor.x = element_line(color = 'grey', size = 0.1),
        panel.border = element_blank(), 
        axis.text.y = element_blank())
f1_5

Stats

  • Repeated-measures ANOVA on log-transformed chases
## 
## Error: as.factor(RealMouseNumber)
##           Df Sum Sq Mean Sq F value Pr(>F)
## Sex        1   0.34  0.3366   0.211  0.648
## Stage      1   0.59  0.5950   0.372  0.544
## Sex:Stage  1   0.01  0.0093   0.006  0.940
## Residuals 80 127.88  1.5984               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value      Pr(>F)    
## Stage      1  5.704   5.704  30.559 0.000000415 ***
## Sex:Stage  1  0.004   0.004   0.023       0.881    
## Residuals 78 14.559   0.187                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Figure 2 - CMS outcomes

Physiological Outcomes

Bodyweight

Stats:

  • Summary:
  • Shapiro-Wilk test

Normality is not violated in any group

  • Levene’s Test
Df F value Pr(>F)
group 3 2.658661 0.0550054
69 NA NA

Groupwise variances are not heteroscedastic (although a trend exists)

Males

Df F value Pr(>F)
group 1 3.419868 0.0742915
30 NA NA

Females

Df F value Pr(>F)
group 1 2.793313 0.1026649
39 NA NA
  • Test

Pairwise

Coat state

Stats:

  • Summary:
  • Shapiro-Wilk test

Normality is violated in both male groups as well as female controls

  • Levene’s Test
Df F value Pr(>F)
group 3 1.535438 0.2130927
69 NA NA

Groupwise variances are not heteroscedastic

  • Transform data:

No transformation helped for the female controls

  • Test

Condition:

Sex:

Condition per sex

Adrenals

Stats:

  • Summary:
  • Shapiro-Wilk test

Normality is not violated in any group

  • Levene’s Test
Df F value Pr(>F)
group 3 0.4508866 0.7174833
69 NA NA

Groupwise variances are not heteroscedastic

  • Test

  • Test

Without interaction:

Pairwise

Cort

Stats:

  • Summary:
  • Shapiro-Wilk test

Normality is violated in all groups

  • Transform the data didn’t help

  • Levene’s Test

Df F value Pr(>F)
group 3 3.555961 0.0186511
69 NA NA

Groupwise variances are heteroscedastic

Per sex Males - not heteroscedastic

Df F value Pr(>F)
group 1 0.4056222 0.5290303
30 NA NA

Females - not heteroscedastic

Df F value Pr(>F)
group 1 0.2798986 0.5997662
39 NA NA
  • Test

Condition:

Sex (not fair, since variances are unequal):

Condition per sex

Correlation plot

Dimensions

## [1] 73 37

Samples

Figure 3 - PCA

Variance explained

PCs vs sex & condition

Check Normality

Normality is not violated in any groups

Check heteroskedasticity

All variances are homogenous

PC1:

##               Df Sum Sq Mean Sq F value Pr(>F)
## Condition      1    3.3   3.304   0.411  0.523
## Sex            1    0.1   0.140   0.017  0.895
## Condition:Sex  1    8.4   8.401   1.046  0.310
## Residuals     69  554.3   8.033

PC2:

##               Df Sum Sq Mean Sq F value         Pr(>F)    
## Condition      1   3.34    3.34   1.233          0.271    
## Sex            1 139.50  139.50  51.449 0.000000000646 ***
## Condition:Sex  1   6.04    6.04   2.226          0.140    
## Residuals     69 187.09    2.71                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PC3:

##               Df Sum Sq Mean Sq F value  Pr(>F)   
## Condition      1  28.70  28.699   7.897 0.00644 **
## Sex            1   1.10   1.097   0.302 0.58449   
## Condition:Sex  1   0.87   0.874   0.241 0.62530   
## Residuals     69 250.74   3.634                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PC4:

##               Df Sum Sq Mean Sq F value  Pr(>F)   
## Condition      1   0.42   0.417   0.159 0.69176   
## Sex            1  27.24  27.237  10.344 0.00198 **
## Condition:Sex  1  30.08  30.076  11.423 0.00120 **
## Residuals     69 181.68   2.633                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PC1

Stats

## 
## Call:
## lm(formula = PC1 ~ DS * Sex, data = pp[pp$Condition == "CMS", 
##     ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.9194 -1.7251 -0.1317  1.6727  5.5421 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   -4.421      2.045  -2.162  0.03754 * 
## DS             2.856      1.287   2.220  0.03301 * 
## SexMale        9.246      2.591   3.568  0.00107 **
## DS:SexMale    -5.834      1.637  -3.563  0.00108 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.554 on 35 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.2857, Adjusted R-squared:  0.2245 
## F-statistic: 4.667 on 3 and 35 DF,  p-value: 0.007594
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## DS           1   6.22    6.22   0.954 0.33533   
## Sex          1   2.30    2.30   0.353 0.55639   
## DS:Sex       1  82.78   82.78  12.695 0.00108 **
## Residuals   35 228.22    6.52                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 9 observations deleted due to missingness

Correlations with behaviors

## [1] 16

Selected features

dat = obj$tests_combined_adj
idVars = c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber', 
           'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props = colnames(dat)[!colnames(dat) %in% idVars]
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                             idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(as.character(MouseID))) %>% 
  rename(Dominance = DS)
d = left_join(dat, baseline)

idVars = c(idVars, 'Dominance')
plots = d  %>% 
  dplyr::select(idVars, props) %>% 
  gather(prop, value, -idVars) %>% 
  group_by(prop) 
plots$prop = plyr::mapvalues(plots$prop, props, props_pretty)
plots$Sex = factor(plots$Sex, levels = c('Female', 'Male'))
plots$Condition = factor(plots$Condition, levels = c('Control', 'CMS'))
plots$SexCondition = paste0(plots$Sex, '_', dat$Condition)
plots$SexCondition = factor(plots$SexCondition, 
                            levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))

selected_props = c('OFT | Distance', 'EPM | Closed arm time')
# selected_props = c('OFT | Inner zone distance', 'EPM | Open arm distance')

res = lapply(selected_props, function(i) {
  pp = subset(plots, prop == i)
  l = c(min(pp$value), max(pp$value))
  p1 = ggplot(pp, aes(x = Condition, y = value, fill = SexCondition)) +
    geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = SexCondition)) +
    geom_point(pch = 21, size = 2, 
               position = position_jitterdodge(jitter.width = 0, 
                                               jitter.height = 0, dodge.width = 0.5)) + 
    scale_fill_manual(values = obj$four_color) + 
    scale_color_manual(values = obj$four_color) + 
    facet_grid(Sex~.) + labs(x = '', y = i) +
    scale_y_continuous(limits = l) +
    theme(legend.position = 'none') + 
    theme(strip.text.y = element_blank())
  
  p2 = ggplot(pp[pp$Condition == 'CMS',], aes(x = Dominance, y = value, fill = Sex, color = Sex)) +
    geom_smooth(method = 'lm', se = F, lty = 2) +
    geom_point(pch = 21, size = 2, color = 'grey20') +
    scale_fill_manual(values = obj$four_color[3:4]) +
    scale_color_manual(values = obj$four_color[3:4]) +
    scale_y_continuous(limits = l) +
    labs(y = '', x = "Baseline David's Score") +
    facet_grid(Sex~Condition) +
    theme(legend.position = 'none', 
          axis.text.y = element_blank(),
          plot.margin = margin(0,0,0,-10, "pt")
    )
  
  cowplot::plot_grid(plotlist = list(p1, p2), axis = 'tb',
                     align = 'h', rel_widths = c(1,1.5))
})
f3_5 = cowplot::plot_grid(plotlist = res, ncol = 2, align = 'v', rel_widths = c(1,1),
                          labels = c('e', 'f'), label_fontface = 'plain')
f3_5

Figure 4 - Rank x CMS

PC1 versus ranks

idVars = c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber', 
           'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props = colnames(dat)[!colnames(dat) %in% idVars]
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                             idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(as.character(MouseID))) %>% 
  rename(Dominance = DS) %>% group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-Dominance), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )

baseline = left_join(baseline, obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])
p <- obj$analyses$PCA_OutcomeVars_BatchAdj %>%
  dplyr::select(RealMouseNumber, PC1)

d = left_join(baseline, p)
d$Sex = factor(d$Sex, levels = c('Female', 'Male'))
d$Condition = factor(d$Condition, levels = c('Control', 'CMS'))
d$SexCondition = paste0(d$Sex, '_', dat$Condition)
d$SexCondition = factor(d$SexCondition, 
                        levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))
d$Rank_3 = plyr::mapvalues(d$Rank, c('alpha', 'beta', 'gamma', 'delta'), 
                           c('alpha', 'beta / gamma', 'beta / gamma', 'delta'))
d$Rank_3 = factor(d$Rank_3, levels = c('alpha', 'beta / gamma', 'delta'))
d$Rank_3_ordinal = as.numeric(as.character(plyr::mapvalues(d$Rank_3,  
                                                           c('alpha', 'beta / gamma', 'delta'), 
                                                           c(1, 0, -1))))
pp = d
pp$value = pp$PC1
l = c(min(pp$value), max(pp$value))

xx = data.frame(
  Sex = c('Female', 'Female', 'Male', 'Male'), 
  x = c(1.5, 2.5, 1.5, 2.5),
  y = c(5, 0, 4.5, 5),
  Sig =  c('*', '', '#', '*'), 
  Line =  c('_', '', '_', '_')
  )
f4_1 = ggplot(pp[pp$Condition == 'CMS',], aes(x = Rank_3, y = value, fill = Rank_3)) +
  geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = Rank_3)) +
  geom_point(pch = 21, size = 2, 
             position = position_jitterdodge(jitter.width = 0, 
                                             jitter.height = 0, dodge.width = 0.5)) + 
  geom_text(data = xx[-3,], inherit.aes = F, aes(x = x, y = y, label = Sig), fontface = 'bold', size = 4,
            color = 'grey30') +
  geom_text(data = xx[3,], inherit.aes = F, aes(x = x, y = y, label = Sig), size = 3, fontface = 'bold', color = 'grey30') +
  
  geom_text(data = xx, inherit.aes = F, aes(x = x, y = y, label = Line), fontface = 'bold', color = 'grey30') +
  
  scale_color_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4])) +
  scale_fill_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4]), 
                    labels = c('\u03B1', paste('\u03B2', '\u03B3', sep = '/'), '\u03B4')) +
  scale_x_discrete(labels = c('\u03B1', paste('\u03B2', '\u03B3', sep = '/'), '\u03B4')) +
  scale_y_continuous(limits = l) +
  facet_grid(Sex~Condition) + labs(x = '', y = 'PC1 score') +
  theme(legend.position = 'none',
        axis.title.x = element_text(family="Helvetica")
  ) 
f4_1

Stats

  • Summary:
  • Shapiro-Wilk test

Normality is not violated

  • Levene’s Test
Df F value Pr(>F)
group 5 0.8531606 0.5225123
33 NA NA

Groupwise variances are homogenous

Test:

Males - omnibus - rank ordinal

Alpha to beta/gamma

Delta to beta/gamma

Females - omnibus

Alpha to beta/gamma

Delta to beta/gamma

Current

Select top PC1 behaviors

  • All features with with a significant adjusted PC1 correlation

All Controls

All CMS

Male Controls

Male CMS

Female Control

Female CMS

Selected props

dat = obj$tests_combined_adj
idVars = c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber', 
           'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props = colnames(dat)[!colnames(dat) %in% idVars]
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                             idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(as.character(MouseID))) %>% 
  rename(Dominance = DS) %>% group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-Dominance), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )

d = left_join(dat, baseline)
idVars = c(idVars, 'Dominance', 'Rank')
d = d  %>% 
  dplyr::select(idVars, props) %>% 
  gather(prop, value, -idVars) %>% 
  group_by(prop) 
d$prop = plyr::mapvalues(d$prop, props, props_pretty)
d$Sex = factor(d$Sex, levels = c('Female', 'Male'))
d$Condition = factor(d$Condition, levels = c('Control', 'CMS'))
d$SexCondition = paste0(d$Sex, '_', dat$Condition)
d$SexCondition = factor(d$SexCondition, 
                        levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))

selected_props = as.character(cors$feature_pretty[cors$Sig == T])
d = d %>% filter(prop %in% selected_props)
d$Rank_3 = plyr::mapvalues(d$Rank, c('alpha', 'beta', 'gamma', 'delta'), 
                           c('alpha', 'beta / gamma', 'beta / gamma', 'delta'))
d$Rank_3 = factor(d$Rank_3, levels = c('alpha', 'beta / gamma', 'delta'))
d$Rank_3_ordinal = as.numeric(as.character(plyr::mapvalues(d$Rank_3,  
                                                           c('alpha', 'beta / gamma', 'delta'), 
                                                           c(1, 0, -1))))
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "

rr = r %>% 
  filter(Sex == 'Male' & comparison == 'alpha_to_control') %>%
  arrange(-p.value) %>% 
  dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)

rr1 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
  geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
  geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
  scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
  coord_flip() + theme(legend.position = 'none') + 
  ggtitle('CMS \u03B1-males vs. Controls')

rr = r %>% 
  filter(Sex == 'Male' & comparison == 'delta_to_control') %>%
  arrange(-p.value) %>% 
  dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
rr2 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
  geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
  geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
  scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
  coord_flip() + theme(legend.position = 'none') + 
  ggtitle('CMS \u03B4-males vs. Controls')

rr = r %>% 
  filter(Sex == 'Female' & comparison == 'alpha_to_control') %>%
  arrange(-p.value) %>% 
  dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
rr3 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
  geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
  geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
  scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
  coord_flip() + theme(legend.position = 'none') +
  ggtitle('CMS \u03B1-females vs. Controls')

cowplot::plot_grid(plotlist = list(rr1, rr2, rr3), ncol = 3)

extra:

Selected prop plots

rr = r %>% filter(prop %in% c('OFT | Distance', 'EPM | Closed arm time')) %>%
# rr = r %>% filter(prop %in% c('EPM | Open arm distance', 'EPM | Closed arm time')) %>% 
  dplyr::mutate(Sig = p.value <=.05)
rr$Sig = ifelse(rr$p.value >.1, NA, 
                 ifelse(rr$p.value >= 0.05, '#', 
                        ifelse(rr$p.value >= .01, '*', 
                               ifelse(rr$p.value >= 0.001, "**", '***'))))

res = d %>% group_by(prop) %>% nest() %>% 
  filter(prop %in% c('OFT | Distance', 'EPM | Closed arm time')) %>%   
  mutate(
    plot = map2(data, prop, function(pp, i) {
      l = c(min(pp$value), max(pp$value))
      p1 = ggplot(pp[pp$Condition == 'Control',], aes(x = as.factor(1), y = value, fill = Rank_3)) +
        geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, fill = 'grey', color = 'grey') +
        geom_point(pch = 21, size = 2, position = position_jitter(height = 0, width = 0)) +
        scale_y_continuous(limits = l) +
        scale_color_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4])) +
        scale_fill_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4])) +
        facet_grid(Sex~Condition) + labs(x = '', y = i) +
        theme(legend.position = 'none',
              strip.text.y = element_blank(),
              axis.text.x = element_blank() 
        ) 
      
      xx = subset(rr, prop == i)
      xx$y = l[2] 
      xx$Rank_3 = plyr::mapvalues(xx$comparison, c('alpha_to_control', 'delta_to_control'), 
                                  c('alpha', 'delta'))
      p2 = ggplot(pp[pp$Condition == 'CMS',], aes(x = Rank_3, y = value, fill = Rank_3)) +
        geom_boxplot(width = 0.5, alpha = 0.5, outlier.shape = NA, aes(color = Rank_3)) +
        geom_point(pch = 21, size = 2, 
                   position = position_jitterdodge(jitter.width = 0, 
                                                   jitter.height = 0, dodge.width = 0.5)) + 
        scale_color_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4])) +
        scale_fill_manual(values = c(obj$colors_rank[1], obj$colors_rank[2], obj$colors_rank[4]), 
                          labels = c('\u03B1', paste('\u03B2', '\u03B3', sep = '/'), '\u03B4')) +
        scale_x_discrete(labels = c('\u03B1', paste('\u03B2', '\u03B3', sep = '/'), '\u03B4')) +
        scale_y_continuous(limits = l) +
        
        geom_text(data = xx, inherit.aes = F, aes(x = Rank_3, y = y, label = Sig), fontface = 'bold', color = 'grey30') +
        
        facet_grid(Sex~Condition) + labs(x = '', y = '') +
        theme(legend.position = 'none',
              plot.margin = margin(0,0,0,-18, "pt"),
              axis.text.y = element_blank(),
              axis.title.x = element_text(family="Helvetica")
        ) 
      cowplot::plot_grid(plotlist = list(p1, p2), axis = 'tb',
                         align = 'h', rel_widths = c(1, 1.3))
    })
  )

Supplement

DS vs SB behaviors

  • Cumulative David’s score based on chases
  • Mean of SB variables over days 2,3,4
wins = 'Chases'
h = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = wins, nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = F)
h = h %>% dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(as.character(MouseID)))
h = left_join(h, obj$mice[,c('RealMouseNumber', 'Sex', 'GroupNumber', 'MouseID')])

p = subset(obj$profiles, !Suspect & Day %in% 2:4) %>% 
  dplyr::select(RealMouseNumber, all_props) %>% 
  group_by(RealMouseNumber) %>% 
  dplyr::summarise_at(.vars = all_props, mean, na.rm = T)

d = left_join(h, p)

idVars = c('RealMouseNumber', 'GroupNumber', 'Sex', 'DS')
props = all_props

p <- d  %>% 
  dplyr::select(idVars, props) %>% 
  gather(prop, value, -idVars) %>% 
  group_by(prop, Sex) %>% nest() %>% 
  mutate(
    M = map(data, function(d) {
      tidy(cor.test(d$value, d$DS, method = 'spearman'))
    })
  ) %>% dplyr::select(-data) %>% unnest() %>% 
  mutate(
    log10pval = -log10(p.value),
    prop = fct_reorder(prop, abs(estimate)), 
    sign = ifelse(estimate > 0, '+', '-'), 
    sig_nominal = p.value <= 0.05
  ) %>% group_by(Sex) %>% 
  mutate(qval = p.adjust(p.value, method = 'BH')) %>% 
  arrange(qval)

plots = d  %>% 
  dplyr::select(idVars, props) %>% 
  gather(prop, value, -idVars) %>% 
  group_by(prop) %>% nest() %>% 
  mutate(
    plot = pmap(list(data, prop), 
                function(d, p) {
                  ggplot(d, aes(x = DS, y = value, fill = Sex)) +
                    geom_smooth(method = 'lm', se = F, lty = 2, aes(col = Sex)) +
                    geom_point(pch = 21, size = 2) +
                    labs(y = p, x = "Baseline David's Score") +
                    facet_grid(.~Sex) +
                    scale_fill_manual(values = obj$colors_sex) +
                    scale_color_manual(values = obj$colors_sex) +
                    theme(legend.position = 'none')
                })) %>% dplyr::select(-data)

p = left_join(p, plots)
  • displaying only behaviors that showed a q-val significant correlation with DS

Hierarchy Panel

  • Hierarchy properties based on chases

  • Daily & mean of days 1:4 and day 5 separately

domProps = c('Steepness', 'Despotism', 'DirectionalConsistency', 'Landau')
domPropsPretty = c('Steepness', 'Despotism', 'Directional Consistency', "Landau's modified h'")

wins = 'Chases'
h = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:5), wins = wins, nPerm = 1,
                      idVars = c('Sex', 'Condition'), verbose = F, daily = T)
h$uGroup = paste0(h$GroupNumber, h$Sex)
plots = lapply(1:(length(domProps)-1), function(i) {
  prop = domProps[i]
  h$prop = h[[prop]]
  lims = c(min(h$prop, na.rm = T), max(h$prop, na.rm = T))
  p1 = ggplot(h, aes(x = as.factor(Day), y = prop, fill = Sex)) +
    geom_line(aes(color = Sex, group = uGroup), lty = 1, size = 0.2, position = position_dodge(0)) +
    geom_point(aes(group = uGroup, col = Sex), pch = 19, size = 1.2, position = position_dodge(0)) +
    labs(x = '', y = domPropsPretty[i])  + 
    scale_color_manual(values = obj$colors_sex) +
    scale_fill_manual(values = obj$colors_sex) +
    scale_y_continuous(limits = lims) +
    theme(legend.position = 'none', axis.title.x=element_blank())
  
  p2 = h %>% dplyr::filter(Day %in% 1:4) %>% 
    group_by(GroupNumber, Condition, Sex) %>% 
    dplyr::summarise(prop = mean(prop, na.rm = T)) %>% 
    ggplot(aes(x = Sex, y = prop, fill = Sex)) +
    geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
    geom_point(pch = 21, size = 2.5, position = position_dodge(0.5)) +
    labs(x = '', y = '')  + 
    scale_y_continuous(limits = lims) +
    scale_fill_manual(values = obj$colors_sex) +
    scale_color_manual(values = obj$colors_sex) +
    theme(legend.position = 'none', 
          axis.text.y = element_blank(),
          axis.title=element_blank()
    )  +
    ggtitle('Baseline')
  
  p3 = h %>% dplyr::filter(Day %in% 5) %>% 
    group_by(GroupNumber, Condition, Sex) %>% 
    dplyr::summarise(prop = mean(prop, na.rm = T)) %>% 
    ggplot(aes(x = Sex, y = prop, fill = Sex)) +
    geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
    geom_point(pch = 21, size = 2.5, position = position_dodge(0.5)) +
    labs(x = '', y = '')  + 
    scale_y_continuous(limits = lims) +
    scale_fill_manual(values = obj$colors_sex) +
    scale_color_manual(values = obj$colors_sex) +
    theme(legend.position = 'none', 
          axis.text.y = element_blank(),
          axis.title=element_blank()
    )  +
    ggtitle('Acute stress')
  cowplot::plot_grid(plotlist = list(p1, p2, p3), ncol = 3, align = 'h',
                     label_fontface = 'plain',
                     rel_widths = c(3.4, 1, 1))
})

lastPlot = lapply(length(domProps), function(i) {
  prop = domProps[i]
  h$prop = h[[prop]]
  lims = c(min(h$prop, na.rm = T), max(h$prop, na.rm = T))
  p1 = ggplot(h, aes(x = as.factor(Day), y = prop, fill = Sex)) +
    geom_boxplot(width = 0.7, alpha = 0.8, outlier.shape = NA, aes(color = Sex)) +
    geom_dotplot(binaxis='y', stackdir='center', dotsize = 1, position = position_dodge(0.7)) +
    labs(x = 'Day', y = domPropsPretty[i])  + 
    scale_fill_manual(values = obj$colors_sex) +
    scale_color_manual(values = obj$colors_sex) +
    scale_y_continuous(limits = lims) +
    theme(legend.position = 'none')
  
  p2 = h %>% dplyr::filter(Day %in% 1:4) %>% 
    group_by(GroupNumber, Condition, Sex) %>% 
    dplyr::summarise(prop = mean(prop, na.rm = T)) %>% 
    ggplot(aes(x = Sex, y = prop, fill = Sex)) +
    geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
    # geom_point(pch = 21, size = 2.5, position = position_dodge(0.5)) +
    geom_dotplot(binaxis='y', stackdir='center', dotsize = 1.3) +
    labs(x = '', y = '')  + 
    scale_y_continuous(limits = lims) +
    scale_fill_manual(values = obj$colors_sex) +
    scale_color_manual(values = obj$colors_sex) +
    theme(legend.position = 'none', axis.text.y = element_blank(), 
          axis.title.y=element_blank()) +
    ggtitle('Baseline')
  
  p3 = h %>% dplyr::filter(Day %in% 5) %>% 
    group_by(GroupNumber, Condition, Sex) %>% 
    dplyr::summarise(prop = mean(prop, na.rm = T)) %>% 
    ggplot(aes(x = Sex, y = prop, fill = Sex)) +
    geom_boxplot(width = 0.5, alpha = 0.8, outlier.shape = NA, aes(col = Sex)) +
    geom_dotplot(binaxis='y', stackdir='center', dotsize = 1.3) +
    labs(x = '', y = '')  + 
    scale_y_continuous(limits = lims) +
    scale_fill_manual(values = obj$colors_sex) +
    scale_color_manual(values = obj$colors_sex) +
    theme(legend.position = 'none', axis.text.y = element_blank(), 
          axis.title.y=element_blank()) + 
    ggtitle('Acute stress')
  
  cowplot::plot_grid(plotlist = list(p1, p2, p3), ncol = 3, align = 'h', 
                     rel_widths = c(3.4, 1, 1))
})
plots[[length(domProps)]] = lastPlot[[1]]
cowplot::plot_grid(plotlist = plots, ncol = 1, align = 'v', axis = 'tblr', 
                   labels = c('a', 'b', 'c', 'd'), label_fontface = 'plain',
                   rel_heights = c(rep(1, length(domProps)-1), 1.1))

Stats:

Normality violated only for Landau’s h in female groups during acute stress

Check heteroskedasticity

Variances are heteroskedastic for despotism

Test

  • Steepness
## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## Sex        1 0.4634  0.4634  21.181 0.000194 ***
## Stage      1 0.0247  0.0247   1.131 0.300839    
## Residuals 19 0.4157  0.0219                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Stage      1 0.00001 0.000015   0.002  0.964
## Sex:Stage  1 0.00109 0.001091   0.157  0.696
## Residuals 19 0.13186 0.006940
  • Despotism
## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value Pr(>F)  
## Sex        1 0.3859  0.3859   7.834 0.0111 *
## Residuals 20 0.9852  0.0493                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Stage      1 0.00707 0.007074   0.757  0.395
## Sex:Stage  1 0.00600 0.006003   0.642  0.432
## Residuals 20 0.18691 0.009345
  • DirectionalConsistency
## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## Sex        1 0.5689  0.5689  24.405 0.0000909 ***
## Stage      1 0.1325  0.1325   5.685    0.0277 *  
## Residuals 19 0.4429  0.0233                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Stage      1 0.00723 0.00723   0.611  0.444
## Sex:Stage  1 0.00270 0.00270   0.228  0.638
## Residuals 19 0.22495 0.01184

Interesting effect of stage

Males:

## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value Pr(>F)
## Stage      1 0.1325 0.13252   3.153  0.114
## Residuals  8 0.3362 0.04202               
## 
## Error: Within
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Stage      1 0.00901 0.009015   0.523   0.49
## Residuals  8 0.13779 0.017223

Females:

## 
## Error: uGroup
##           Df Sum Sq  Mean Sq F value Pr(>F)
## Residuals 11 0.1067 0.009702               
## 
## Error: Within
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Stage      1 0.00092 0.000916   0.116   0.74
## Residuals 11 0.08716 0.007924

  • Landau
## 
## Error: uGroup
##           Df Sum Sq Mean Sq F value Pr(>F)
## Sex        1 0.1363 0.13632   2.113  0.162
## Stage      1 0.0709 0.07092   1.099  0.308
## Residuals 19 1.2260 0.06452               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Stage      1 0.0012 0.00121   0.050  0.825
## Sex:Stage  1 0.0384 0.03841   1.598  0.222
## Residuals 19 0.4568 0.02404

Baseline non-parametric:

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Landau by Sex
## Kruskal-Wallis chi-squared = 5.8385, df = 1, p-value = 0.01568

Acute stress non-parametric:

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Landau by Sex
## Kruskal-Wallis chi-squared = 0.022159, df = 1, p-value = 0.8817

Alpha/Delta-males and Alpha females

dat = obj$tests_combined_adj
idVars = c('Color', 'Sex', 'Condition', 'Batch', 'RealMouseNumber', 
           'GroupNumber', 'MouseID', 'GroupID', 'uGroup', 'GroupType')
props = colnames(dat)[!colnames(dat) %in% idVars]
baseline = hierarchyProfiler(subset(obj$master, !Suspect & Day %in% 1:4), wins = 'Chases', nPerm = 1,
                             idVars = c('Sex', 'Condition'), verbose = F, daily = F) %>% 
  dplyr::select(GroupNumber, Sex, Condition, DS) %>% unnest() %>% 
  mutate(MouseID = as.numeric(as.character(MouseID))) %>% 
  rename(Dominance = DS) %>% group_by(GroupNumber, Sex) %>% 
  mutate(
    Rank = rank(-Dominance), 
    Rank = plyr::mapvalues(Rank, 1:4, c('alpha', 'beta', 'gamma', 'delta')) ,
    Rank = factor(Rank, levels = c('alpha', 'beta', 'gamma', 'delta'))
  )

d = left_join(dat, baseline)
idVars = c(idVars, 'Dominance', 'Rank')
d = d  %>% 
  dplyr::select(idVars, props) %>% 
  gather(prop, value, -idVars) %>% 
  group_by(prop) 
d$prop = plyr::mapvalues(d$prop, props, props_pretty)
d$Sex = factor(d$Sex, levels = c('Female', 'Male'))
d$Condition = factor(d$Condition, levels = c('Control', 'CMS'))
d$SexCondition = paste0(d$Sex, '_', dat$Condition)
d$SexCondition = factor(d$SexCondition, 
                        levels = c('Female_Control', 'Male_Control', 'Female_CMS', 'Male_CMS'))

selected_props = as.character(cors$feature_pretty[cors$Sig == T])
d = d %>% filter(prop %in% selected_props)
d$Rank_3 = plyr::mapvalues(d$Rank, c('alpha', 'beta', 'gamma', 'delta'), 
                           c('alpha', 'beta / gamma', 'beta / gamma', 'delta'))
d$Rank_3 = factor(d$Rank_3, levels = c('alpha', 'beta / gamma', 'delta'))
d$Rank_3_ordinal = as.numeric(as.character(plyr::mapvalues(d$Rank_3,  
                                                           c('alpha', 'beta / gamma', 'delta'), 
                                                           c(1, 0, -1))))

r1 = d %>% group_by(prop, Sex) %>% 
  nest() %>% 
  mutate(
    res = map(data, function(dds) {
      dds = dds[!(dds$Condition == 'CMS' & dds$Rank_3_ordinal %in% c(0, -1)),]
      
      # dds = dds %>%
      #   ungroup() %>%
      #   dplyr::mutate(
      #     value_ranked = rank(value)
      #   )
      # # x = as.data.frame(summary(lmPerm::aovp(value_ranked ~ Condition, dds))[[1]])
      
      x = as.data.frame(summary(lmPerm::aovp(value ~ Condition, dds))[[1]])
      x = x[1,]
      
      xx = dds %>% group_by(Condition) %>% 
        dplyr::summarise(M = mean(value, na.rm = T))
      x$direction = ifelse(xx$M[xx$Condition == 'CMS'] > xx$M[xx$Condition == 'Control'], 'positive', 'negative')
      # x = tidy(kruskal.test(value ~ Condition, dds))
      x$comparison = 'alpha_to_control'
      return(x)
    })
  ) %>% dplyr::select(-data) %>% unnest()
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
## [1] "Settings:  unique SS "
rr = r %>% 
  filter(Sex == 'Male' & comparison == 'alpha_to_control') %>%
  arrange(-p.value) %>% 
  dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)

rr1 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
  geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
  geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
  scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
  coord_flip() + theme(legend.position = 'none') + 
  ggtitle('CMS \u03B1-males vs. Controls')

rr = r %>% 
  filter(Sex == 'Male' & comparison == 'delta_to_control') %>%
  arrange(-p.value) %>% 
  dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
rr2 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
  geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
  geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
  scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
  coord_flip() + theme(legend.position = 'none') + 
  ggtitle('CMS \u03B4-males vs. Controls')

rr = r %>% 
  filter(Sex == 'Female' & comparison == 'alpha_to_control') %>%
  arrange(-p.value) %>% 
  dplyr::mutate(sig = p.value <= 0.05)
rr$prop = factor(rr$prop, levels = rr$prop)
rr3 = ggplot(rr, aes(x = prop, y = -log10(p.value))) +
  geom_hline(yintercept = -log10(.05), color = 'grey', lty = 2) +
  geom_point(aes(fill = sig), pch = 21, color = 'black') + labs(y = '-log10(p-value)', x = '') +
  scale_fill_manual(values = c('grey', hm.palette(9)[1]), name = 'Significant') +
  coord_flip() + theme(legend.position = 'none') +
  ggtitle('CMS \u03B1-females vs. Controls')

cowplot::plot_grid(plotlist = list(rr1, rr2, rr3), ncol = 3)